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Abstract 

Three coarse-grained models of the double-stranded DNA are pro- 
posed and compared in the context of mechanical manipulation such 
as twisting and various schemes of stretching. The models differ in the 
number of effective beads (between two and five) representing each nu- 
cleotide. They all show similar behavior and, in particular, lead to a 
torque-force phase diagrams qualitatively consistent with experiments 
and all-atom simulations. 
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1 Introduction 



Manipulation of large biomolecules by means of atomic force mi- 
croscopy, optical tweezers, and other nanotechnological devices is play- 
ing an increasingly growing role in elucidating mechanisms of bio- 
logically relevant processes [U El E]. The dynamical data obtained 
through mechanical manipulation usually requires theoretical interpre- 
tation that can be reached through numerical simulations. This need is 
especially apparent when dealing with proteins (see, e.g., [4]) because 
of their strongly inhomogeneous network of interactions between the 
amino acids. The inhomogeneity results, for instance, in a protein- 
dependent pattern of peaks when the force of resistance to pulling at 
constant speed is plotted against elongation (see, e.g. [5j [6j [7j [8]). 

All-atom simulations have contributed to understanding of the 
large conformational changes in proteins induced by the mechanical 
manipulation (see e.g. [9j \W\ CD]). However, such simulations are 
inherently limited by the time scales and system sizes that can be 
studied. One way out is to accelerate the modeled processes by orders 
of magnitude relative to their experimental realizations which may, 
undesirably, distort the physics involved. Another way is to use coarse 
grained models. These models can also be applied to other processes 
involving large conformational changes such as folding or thermal un- 
folding of biomolecules. They may also find applications in studies of 
systems which are much larger in size. 

In this paper, we focus on coarse grained models describing the 
double stranded DNA (dsDNA). The simplest coarse graining scheme 
involves treating dsDNA as a polymer endowed with a local stiffness 
|12j . The resulting model seems to be appropriate in studies of stretch- 
ing at high temperatures, when the entropic effects dominate, and in 
studies of phenomena related to fluid flow and the hydrodynamic in- 
teractions [13]. The corresponding characteristics dimension in this 
model can be measured by the Kuhn length or the hydrodynamic ra- 
dius. The typical values of these parameters are of order 106 and 77 
nm [13] respectively which encompasses more than 200 distances (of 
0.4 nm ) between successive pairs of nucleotides in the dsDNA. This 
level of coarse graining is thus too crude to be useful when describing 
mechanical manipulation of the system where nanometric features are 
detectable. Here, we discuss coarse grained models with structures 
that are resolved at the level of a single nucleotide. 

There are a variety of possible ways to manipulate the DNA duplex. 
The ones that we consider in this paper are schematically represented 
in Figure 1. Unzipping, whether at constant speed or at constant 
force, corresponds to scheme A in this Figure. The unzipping process 
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involves breaking of one hydrogen bond at a time resulting in a force, 
F, of order 13-15 pN, that rapidly ondulates with the amplitude of 
order 1 pN as the pulling distance, d, is increased [HI EES]. The force 
pattern has also been found to depend on the pulling speed, v p , weakly 
|14j . but would be expected to depend on temperature, T, more sub- 
stantially. Even though the variations in the F — d patterns when 
pulling at constant speed take place on the nucleotidic length scales, 
their interpretation in terms of the specifics of the sequence is difficult 
because of the noise due to thermal fluctuations. Recently, Baldazzi 
et. al. [16] has suggested, however, that if instead the mechanical un- 
zipping is performed at constant force then Bayesian methods of the 
corresponding sequence prediction should be nearly error-free. Meth- 
ods for extracting kinetic information from constant force experiments 
have been discussed recently in the context of DNA unzipping in a 
nanopore |17j . 

Schemes B and C involve stretching at the opposite ends of the 
duplex. The two schemes employ distinct mechanisms of resistance 
to stretching. Mechanism B involves shear which is responsible for 
generation of the strongest force clamps in biomolecules [4] and the 
maximum force obtained depends on the number of the bonds that 
are sheared simultaneously. Mechanism C, on the other hand, leads 
to localized unravelling and generates a size-independent force. When 
using micropipettes, as in the experiment by Cluzel et al. [18J . one 
probably combines schemes B and C. The resulting F — d pattern has 
three stages: one starts off with a long period of a nearly constant 
force, after which there is a steady increase (in the 120 pN range) 
which finally is followed by a sudden drop to zero. A similar pattern 
of behavior arises in simulations involving anisotropic pressure |19| 
(these studies have been performed for a dsDNA with about 10 pairs 
of nucleotides) . 

Still another way of manipulating the dsDNA has been employed by 
Oroszi et al. [20] and it involves applying a torque, G. Wereszczynski 
and Andricioaei [21] have generalized it still further, in their all-atom 
simulations, by considering a simultaneous application of a force and 
torque, as shown in scheme D in Figure 1. They have predicted ex- 
istence of a rich phase diagram of possible structures on the F — G 
plane. There have been experimental studies involving torque pro- 
duced in an optical trap [22j [20j [23]. Bryant et. al. [22] have found 
that the torque needed to transform the B-DNA conformation into 
the left-hand twisted L-DNA form is G = —9.6 pN nm (the negative 
sign indicates twisting against the native sense of turn in the dsDNA), 
while to transform it further into the Pauling-like (P-DNA) form the 
required torque is 34 pN nm. They have also observed that by pulling 
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the dsDNA molecule with the force of about 65 pN one induces a tran- 
sition from the B-DNA into the S-DNA form which is 60% longer then 
the B-DNA. 

In this paper, we construct three variants of the nucleotide-based 
coarse grained dynamical models of the DNA duplex. The coarse grain- 
ing method introduces several effective objects, referred here as beads, 
that represent a single nucleotide. The models differ primarily by the 
number of beads involved. We compare the workings of the three 
models for a 22-base-pair system, or shorter, and use them to eluci- 
date the mechanisms of rupture in processes corresponding to schemes 
A through D. One conclusion of our studies is that even though var- 
ious dynamical details differ between the models, all of them can be 
considered adequate and ready to be applied to larger systems. In par- 
ticular, all of the coarse grained models studied lead to a transitions 
of the usual right-hand-twisted B-dsDNA form to the L-dsDNA form 
and to the P-DNA form on application of a appropriate torque. 

2 An overview of the models used 

Physical properties of the DNA double helix are quite distinct [24] 
compared to other biomolecules. Its strong stiffness comes from the 
braided nature of its structure combined with the presence of the base 
stacking interactions. Furthermore, the phosphate groups in the DNA 
backbone carry substantial electric charges. All of these features are 
employed by the cell's machinery in the processes of copying, transcrib- 
ing and packaging of the DNA. For example, helicases which unwind 
the double helix to provide single-stranded templates for polymerases, 
have evolved as motors that are capable of moving along the torsion- 
ally constrained DNA molecules. Topoisomerases break and reconnect 
the DNA to relieve a torsional strain that accumulates ahead of the 
replication fork. Finally, the DNA-binding proteins get docked to the 
DNA by means of guidance mechanisms which seem to be primarily 
electrostatic in nature. 

Our models address the mechanical properties of the dsDNA and 
do not aim at determining the electrostatic potential outside of the 
duplex. The models are built in analogy to the Go-like models of 
proteins, especially in the specific implementation proposed in refs.[25l 
EZJEElll]- I n the case of proteins, the model represents the system by 
its C a atoms which are tethered together by harmonic interactions. 
The native contacts, such as the hydrogen bonds, are described by 
the Lennard-Jones potentials. The Langevin overdamped thermostat 
with random forces mimics fluctuational effects of the solvent. 61 other 
variants of this basic model of a protein are discussed and compared 
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in ref. [29]. 

The dsDNA has a simpler elastic structure than proteins since a 
pair of nucleotides can bind only in two ways: either by forming two 
(A-T) or three (G-C) hydrogen bonds. When trying to build a coarse 
grained model for a DNA one is first inclined to assign a single bead to 
a nucleotide and to locate it at the phosphorus (P) atom. This may be 
acceptable for a single strand DNA, provided the local chain stiffness 
terms are included. However, for the dsDNA this procedure would 
lead to a distance of 17 A between the P atoms in a hydrogen bonded 
pair of the nucleotides. Such a relatively large distance would intro- 
duce too much mechanical instability in the model but appears to be 
adequate to study conformational changes in dsDNA nanocircles and 
submicron-sized plasmids with torsional stress [30]. A more detailed 
approach, denoted here as model I, involves representing the A and T 
nucleotides by four beads and the G and C nucleotides by five beads. 
One of the beads represents the phosphate group, another the sugar 
group, and the remaining beads participate in formation of either two 
or three hydrogen bonds, depending on the specificity. The hydrogen 
bond interactions are represented by the effective Lennard- Jones po- 
tentials and other bonds, being structural in nature, are described by 
the harmonic potentials with large elastic constants. The schematic 
construction of this model is shown in Figure 2. 

More simplified approaches involve reduction in the number of 
beads representing each nucleotide. In a model denoted here as model 
II, we mimic the ribose-phosphate groups by one bead and the base 
by another bead as illustrated by the top of Figure 3. In this model, 
the distinction between the A:T and G:C pairing interactions comes 
not through introduction of separate contacts for each hydrogen bond 
but through adjustment of the amplitude of the effective base-base 
Lennard- Jones potential by a factor of 2 or 3 respectively. 

In between models I and II there is another model, denoted here 
as model III, that has been introduced by Knott et al. [31] in the 
context of the salt induced melting. Model III involves three beads 
as illustrated in the bottom part of Figure 3. The beads represent 
the phosphate, sugar, and base groups correspondingly. In model III, 
the backbone chain of one strand is constructed by linking the sugar 
group to the P atom on the same nucleotide and to the P atom on a 
preceding nucleotide. Thus the backbone chain has the appearance of 
a zigzag line. In contrast, in models I and II, the backbone chain is 
formed by tethering the consecutive P atoms. 

As a model system we consider the structure coded as 
119D in the Protein Data Bank which has been determined by 
Leonard and Hunter [32]. It corresponds to the sequence 5'- 
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D(CGTAGATCTACGTAGATCTACG)-3'. The ground state confor- 
mations of this system as captured by the three models are shown in 
Figure HI Longer structures can be obtained, for instance, by repeat- 
ing this basic unit, or by constructing a synthetic sequence-dependent 
structure that makes use of average geometrical parameters associated 
with a single base pair. 

3 Model I: the 4- or 5-bead description 

We start by introducing three different types of beads, p, h, and 
b as illustrated in Figure [2l The p-beads are meant to represent the 
backbone which is made of the phosphate groups. The p-bead is placed 
at the position of the C4* atom in the molecule of ribose. This place- 
ment achieves two goals. First, it represents the DNA phosphate chain 
by the p-beads. Second, it locates the p-bead close to the base beads. 
The C4* atom is the ribose ring atom that is closest to the phosphate 
group. The h-beads represent the 'head' atoms which may act either 
as donors or as acceptors in the hydrogen bonds. In the C nucleotides, 
the h-beads are located on the 02, N3, and N4 atoms of the bases. In 
the G nucleotides, - on the 06, Nl, and N2 atoms. Finally, in the A 
and T nucleotides - on the N6, Nl, and N3, 04 atoms respectively. 
The h-beads are linked to their 'bases', i.e. to the supporting b-beads. 
In the native state, the b-beads are located half-way between the p- 
bead and the center of mass of the h-beads at each nucleotide. The 
overall scheme results in about 4- to 5-fold reduction in the number of 
the degrees of freedoms compared to the all-atom approach. 

In Table Q] we provide the Cartesian coordinates of the beads that 
define the model. Knowledge of these coordinates should allow for 
a generation of a synthetic dsDNA based on the sequence. The co- 
ordinates have been obtained by making averages of the geometric 
parameters in the PDB structure corresponding to the basic sequence 
studied. 

The p-beads are tethered into two separate chains and thus form 
two backbones. The tethering is accomplished through the elastic 
potential 

V[f +1 = K b ■ (f Pi - r Pi+1 - d pm+1 ) 2 , (1) 

where index i enumerates consecutive nucleotides and d PiPi+1 are the 
distances between the consecutive p-beads in the native state. These 
distances vary from bead to bead. Their mean is equal to 5.8 A and the 
standard deviation is close to 0.3A. The mean geometrical parameters 
cited in this description can be used in a general sequence-dependent 
construction of a synthetic dsDNA structure. The elastic constant is 
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taken to be equal to = 50 eA -2 , where e is the energy scale corre- 
sponding to the internucleotidic hydrogen bonds, as defined below. 

The steric constraints of the DNA sugar-phosphate backbone are 
represented by the following two potentials for the bond and dihedral 
angles of p-beads' backbone. 

M-2 

V B =J2 Ke{0 t -6 m f , (2) 

i=Ni 

AT-3 

V D =J2 + cos (^ ~ M) + Kfcl + cos3(& - fa))} . (3) 

i=l 

The bond angle 0i is measured between the pi — p%+\ and pi + \ — pi+2 
bonds, and the dihedral angle (pi is an angle between two planes: one 
of them is determined by the pi — Pi+i and pi + \ — Pi+2 bonds, and the 
second one by the Pi+\ — Pi+2 and pi + 2 — Pi+3 bonds, the subscript 
indicates the native values, J\f denotes the number of nucleotides in 
one chain. We take Kg = 20e/(rad) 2 , = l.Oe, K'^ = 0.5e in analogy 
to ref. [33]. 

All of the interbead interactions within one nucleotide are taken to 
be harmonic so that the corresponding potentials read 

v if = E K * ■ (*k - ^ - d n»i ) 2 ' ( 4 ) 

U,fJ, 

where the indices v and \i label beads belonging to the i'th nucleotide. 
The equilibrium distances d^ v take values as in the native structure and 
they range from 2.3 ± 0.1 A for the neighboring h-beads (according 
to notation in Figure [2j pairs: h\j-hij and h2j-hsj) to 3.6 ± 0.7 A 
between the b-bead and h-beads. 

The h-beads on one chain are capable of making hydrogen-bond 
contacts with the h-beads on the opposite chain. In the simplest 
version, we follow the prescription used previously for proteins and 
describe these contacts by the Lennard-Jones potential 

V hh = Ae[{ ^ ) 12_ ( ^Hhl ) 6 ] ^ (5) 
Vfiihj Thill j 

where i and j are the paired residues and r^h = l*"^ — r^.|. The 
parameters a^hj are chosen so that each contact in the native con- 
formation is stabilized in the minimum of the potential. Essentially, 
(Jhihj = ^~ l ^dhihj- The value of d^hj, the distance between h-beads 
making bond, is equal to 2.66 ± 0.14 A. For proteins, the choice of 
the form of the contact potential has turned out to be of much less 
importance than the correct determination of the contact map |29j . 
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The stacking interactions between consecutive b-beads in each 
chain are also accounted for by the Lennard-Jones interactions: 



V$ +1 = 4e[(^±l) 12 - (^±i) 6 ] . (6) 

r bib i+1 r bib i+1 

The distance between the stacking pairs of the b-beads is 4.43 ± 0.42 
A in the native structure. 

All of the interactions discussed above arise in the native state. 
However, distorted conformations may lead to new interactions. In 
the spirit of the Go-like models, we describe these by purely repulsive 
potentials (the Lennard-Jones potentials which are cut at the minimum 
and shifted). The hard sphere diameters of the h- b- and p-beads are 
taken to be equal to 2.0 A, 3.4 A, and 6.0 A respectively. The large 
effective size of the p-bead prevents the chains from crossing and self- 
crossing. 

All beads are endowed with the same mass, m, and the equation 
of motion of each is described by the Langevin equation 

mf = -7(f) + F c + f , (7) 

which provides thermostating and mimics dynamical effects of the sol- 
vent. Here r is the position of the bead, F c is the net force on it due to 
potentials, 7 is the friction coefficient, and T is a white noise term with 
the dispersion of \J2'ykBT, where &b is the Boltzmann constant and 
T is the temperature. The dimensionless temperature, ksT/e, will be 
denoted by T. The friction coefficient 7 is equal to 2m/ V where r is a 
characteristic time scale. The dynamics are meant to be overdamped 
so the characteristic time scale corresponds to a diffusional passage 
of a molecular distance (~ 3 A) and is thus of order 1 ns. For small 
damping, r would correspond to time scale of (ballistic) oscillations in 
the Lennard-Jones well which is significantly shorter. The equations 
of motion are solved by the fifth order predictor-corrector scheme [34J. 

As the average value of energy for hydrogen bond interaction in 
dsDNA we chose 0.6 kcal/mole, while in [35] it was chosen around 0.5- 
0.7 kcal/mole, and in [31]: 0.66 kcal/mole. On average, 2.5 hydrogen 
bonds are created between the bases in dsDNA. Hence the total average 
energy of interaction between paired bases in the DNA helix is about 
1.5 kcal/mole in our models. This choice is consistent with T = 0.4 
corresponding to T=300 K. The corresponding unit of the force, e/A 
should be then of order 100 pN. In the entropic limit, the hydrogen 
bond potentials matter much less than the thermal fluctuations. In 
our model, this starts to happen at T of about 0.5 - 0.6. 



4 Model II: the 2-bead model descrip- 
tion 



In the 2-bead model, we consider beads denoted by p and b at 
each nucleotide as shown in Figure El The p-beads are placed in po- 
sitions of the C4* atom and mimic the phosphate-ribose chain of the 
DNA molecule. The harmonic tethering potential, as well as the bond 
and the dihedral angle potentials are introduced in analogy to model 
I. The b-bead in each nucleotide is placed in the geometrical center 
of the base. The absence of the h-beads of model I is compensated 
by introducing hydrogen-bond-like interactions between the b-beads. 
The average Cartesian coordinates of the model beads are provided in 
Table [TJ In order to distinguish between the A:T and G:C pairs in 
the DNA sequence, we strengthen the amplitude of the correspond- 
ing Lennard-Jones potential by the factor of 2 or 3 respectively. The 
stacking potential, between the neighbouring b-beads along each DNA 
is described as in model I. The hard sphere diameters of the beads re- 
main defined as in model I. Average value of the bead mass is m = 162 
g/mole and the distance between the beads that make effective hy- 
drogen bonds is 5.5 ± 0.8. We have deteremined that the persistence 
length in model II at T = 0.4 is about 50 nm. 



5 Model III: the 3-bead description 

Model III, introduced in ref. [3T] and shown schematically in Fig- 
ure provides a 3-bead description. It differs from model II primarily 
as a result of a different treatment of the backbone chain. In model II, 
the sugar and phosphate groups of the same nucleotide are represented 
by one bead, whereas in model III the two groups are represented by 
separate beads so that the backbone chain is formed by connecting 
sugar bead (s) of one nucleotide to the phosphate bead of the nu- 
cleotide that follows in the sequence. 

The distinction between the phosphate and sugar groups is impor- 
tant in this model because it facilitates introduction of electrostatic 
charges on the phosphate beads. The charges are introduced to de- 
scribe interactions of the DNA with ions in the solvent but they also 
affect the p-p distances through the resulting Coulombic repulsion. 
The corresponding potential is given by V e lec,%j = ^l^lr- ■ e~ Vi ^ KD , 
where kd is the Debye constant, and its value depends on the ionic 
strength of the solution. For standard ionic strengths, kd ranges from 
11 to 15 A (eg. when [Na 2+ ] = 50mM, one obtains kd = 13. 6 A). Here, 
we do not take this term into account, as its effect on the p-p distances 
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is minor and because our focus is on mechanical manipulations and not 
on the effects resulting from variations of the ionic strength. 

Other potentials used in this model are analogous to those used 
in models I and II. The exception is the base pairing potential. We 
describe it by the effective Lennard-Jones potential whereas Knotts et 
al. [HI] by the 10-12 potential 

V b P= V 4e fep [5(^) 12 -6(^) 10 ] , (8) 

base pairs J J 

where e^p depends on the type of base pair (AT or GC, while ecc = 
3/2 • ecc), and cr^ is around 2.9 A for all paired bases. 

Hyeon and Thirumalai [35] have recently considered a model of the 
RNA hairpin in which every nucleotide is represented by three beads 
which correspond to the phosphate, sugar and base groups which is 
analogous to model III of the double helix and to the model of Knotts 
et al. [HI]. Similar to ref. [HI], the Debye-Hueckel potential between 
the phosphate beads is introduced to account for screening by con- 
densed counterions and for the hydratation effects. However, there are 
differences pertaining to the nonbonded potentials. In addition to the 
base pairing potential, Hyeon and Thirumalai introduce a possibility 
of stacking interactions between the base beads. Such interactions do 
not arise between the successive nucleotides, but may arise in, e.g., the 
head of an RNA hairpin. These stacking interactions are responsible 
for existence of the more complicated conformations, like the hairpin, 
that the RNA may adopt. As to the base pairing potential, Hyeon 
and Thirumalai take the Lennard-Jones potential without making a 
distinction between the number of the hydrogen bonds involved. It 
affects the base beads which are within the distance of 7 A and the 
corresponding depth of the potential is 1.8 kcal/mole. Non-native base 
bead interactions are repulsive and correspond to the energy parameter 
of 1.0 kcal/mole. 

Another simple model of the DNA has been recently proposed by 
Ouldridge et al. [36] in the context of self-assembly of DNA nanos- 
tructures in which the twisting character of the individual strands is 
disregarded. In this model, every nucleotide is represented by a softly 
repulsive sphere of diameter 1=6.3 A and by another smaller sphere 
attached (nearly) rigidly to the center of the repulsive sphere at a 
distance of 0.3Z away from it and perpendicularly to the backbone. 
The smaller sphere provides a center of attraction to another small 
sphere and thus plays the role of a base in the DNA strand. Four 
types of the base site are considered and the Lennard-Jones attraction 
links the complementary bases. Additionally, the model incorporates 
a monomer-to-monomer bending energy to provide stiffness. In the 
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ground state, two strands run parallel to each other like in a (3 sheet 
in proteins. 

Throughout the paper, we use the open square, solid square, and 
triangle symbols to denote results corresponding to models I, II, and 
III respectively. 

A convenient way to characterize the unravelling process is by pro- 
viding the distance at which a given contact breaks down for the last 
time. A contact is said to be broken if the corresponding distance 
exceeds 1.5 times the length parameter a in the Lennard- Jones poten- 
tial associated with the contact. For models II and III, where there is 
only one connection between paired bases (through the b-beads) the 
contacts are labelled by the nucleotide number, Z, as counted from the 
side that is being pulled. In the case of model I, where each nucleotide 
can participate in either 2 or 3 contacts, the graphical representation 
uses the ordinates of Z — 0.1 and Z + 0.1 for residues A and T whereas 
it involves I — 0.2, Z and Z + 0.2 for residues G and C. 

6 Unzipping at constant speed in 
scheme A 

The stretching scheme A shown in Figure Q] leads to the unzipping 
process in which the hydrogen bonds break by starting from the end 
that is being pulled. These bonds are enumerated by the index Z which 
is being counted from the pulling end. Figure [5] shows the force vs. 
displacement patterns at T=0 for two different values of the pulling 
velocity, v p , of 0.05 and 0.005 A/r and for the three models discussed. 
At temperatures that are lower than 0.12 e/fcg, the F — d patterns 
are qualitatively similar to the ones shown in Figure [5] in the sense 
that the individual force peaks can usually be related to unravelling of 
specific base pairs in the DNA sequence. 

In the initial stages of unravelling, all bonds that exist in the system 
get adjusted to some extent and it is only later on that the unravelling 
process becomes more site specific. This is well seen in model I for 
d < 50 A, especially for the higher pulling speed, where all force peaks 
are quite similar. In the other two models, this transient distance is 
much shorter because the number of the adjustable bonds is smaller. 
However, after this transient stage is passed, one can read off the pair 
sequence of the double helix from the F — d patterns easily because 
the higher peaks arise due to breaking of the G:C bonds. The smaller 
the pulling speed the crisper the recognition of the sequence. 

At T = 0.4, we observe the peaks with the maximal value of 0.39 
e/A (which corresponds to about 40 pN) - Figure [6l These values are 
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about 20 pN higher than the experimental results [HI [15], where for 
v p = 200 nm/s one obtains peaks of 18-20 pN 

The bottom panels of Figures [5] and [6] show the scenario diagrams 
arising in scheme A. At sufficiently low temperatures the zipping 
process is seen to be proceeding linearly in time with minor differ- 
ences in slopes between the models. At higher temperatures, like for 
T = 0.4e/A in model I, the scenario diagram data points acquire a 
curved appearance. This indicates that the thermal fluctuations rup- 
ture bonds at the idle (and not anchored) end of the dsDNA before 
the mechanical unzipping process gets to them. 

For T = the force peaks have values of 1-2 e/A in all three models. 
On increasing the temperature, the peaks and < F >, i.e. the force 
averaged over the duration of the full unravelling process (till F drops 
to 0), get lowered in a monotonic fashion, as shown in Figure [7J which 
presents the results for unravelling with velocities v p = 0.05 A/r and 
v p = 0.005 A/r. At the higher temperatures, model II yields the 
biggest mean forces, independent of the pulling speeds. Models I and 
III yield comparable mean forces at these temperatures and it depends 
on the velocity which one is stronger of the two. 

With lowering the unravelling velocity, the calculated average forces 
decrease. In Figure [8] we present the average forces for model I for 
three different unravelling velocities. There is a big decrease of the 
< F > between the v p = 0.05 A/r and v p = 0.01 A/r, while the 
results obtained for v p = 0.01 A/r and v p = 0.005 A/r are close to 
each other. Generally, the smaller the velocity, the smaller the < F > 
values. The similar dependencies are observed for models II and III. 
Around T = 0.3 e/A the rapid decrease in < F > switches to a nearly 
constant behavior at higher temperatures. 

7 Stretching at constant speed in 
scheme B 

The B-type stretching has an entirely different nature than the A- 
type one. In this scheme, only one chain undergoes active stretching 
and this effect in turn influences the companion chain through the 
hydrogen-bond contacts. FigureEJillustrates the mechanics of this kind 
of manipulation. The snapshots (as obtained within model II) show 
that a full extension of one chain results in a substantial distortion of 
the other chain. Furthermore, F depends on d in a monotonic fashion. 
There are no force peaks even at T=0. The F — d curves for models 
I and II coincide and a bigger force arises faster than in model III 
because of a more direct transmission of tension between the p-beads. 
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The scenario diagrams also look distinct compared to scheme A 
and in a way which is more sensitive to T. At T=0 many contacts 
are broken nearly simultaneously. At finite temperatures, the contacts 
at the extremities get ruptured before unravelling of the contacts in 
the middle in each of the models studied. The higher the T, the 
earlier particular contacts break down. We have observed insignificant 
dependence of the rupture distances on the pulling velocity. The whole 
process results in unravelling both of the hydrogen-bond and of the 
stacking contacts. 

8 Stretching at constant speed in 
scheme C 

In the C-type stretching, one chain is made to slide along its com- 
panion until the two chains separate as shown in Figure \W\ The F — d 
curves display a major peak which is an order of magnitude larger than 
the force peaks observed in scheme A in all of the three models studied. 
The emergence of this major peak is due to an increasingly cooperative 
resistance to manipulation of many contacts that are sheared simul- 
taneously. Once the rupture takes place, the force drops down to the 
level corresponding merely to the thermal noise. The cooperation level 
appears to be the greatest in model III, followed by model II. In each 
of the models, the rupture of all contacts is nearly simultaneous. 

The maximal force peak dependence on temperature shown in Fig- 
ure [TT] for two velocities indicated is similar to what was observed for 
< F > in scheme A, especially at low temperatures. Models II and 
III yields are found to yield comparable forces which are also notice- 
ably larger than in model I except at the low temperature end. For 
the C-type stretching, the dependence on the stretching velocity is 
weak as demonstrated in Figure [T2l The inset demonstrates that the 
dependence is nearly logarithmic. 

9 Stretching at constant force and con- 
stant torque in scheme D 

We now consider the tensile and torsional manipulations of the 
dsDNA and focus on the determination of the corresponding force- 
torque phase diagram. 

We introduce the torsional stress of the dsDNA molecule in the 
following way. At one end of the molecule we choose two vectors 
defining the plane. First vector is defined by the positions of the 
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extreme p-beads at the chosen end of the dsDNA. The second vector 
defining the plane is a cross product of the first vector and the dsDNA 
axis (which is defined by the midpoints of the extreme beads on both 
ends of DNA molecule). In the plane constructed in the way described 
above we add two more beads as shown in Figure H3] so that a square 
frame of four beads is formed. All beads in this frame are connected 
by the springs as in other structural bonds. The extreme beads on the 
other end of the molecule are anchored at their starting positions. 

The torsion is applied to the DNA molecule by application of a 
force to each of the four beads in the square frame. The torque is 
perpendicular to the frame. The torsion is considered positive if it 
agrees with the sense of the twist in the double helix in the B form 
and it is considered negative otherwise. The stretching force F is a 
resultant force applied to all of the four beads in the frame and along 
the dsDNA axis. 

We first consider model II at two temperatures: 0.2 and 0.4 e/A. 
The results are presented in the phase diagram in the Figure [HI The 
boundaries of the phases are approximate and are quite similar in both 
models. These phase diagrams are also similar in appearance to those 
established experimentally [22], [2(2 [23] and theoretically [39l [22] , 

We start from the dsDNA B-form structure and observe the tran- 
sitions into other phases of the DNA structure. At T = 0.4 e/A, the 
DNA transforms into the L-form under the torque G of around —1.5 e. 
At the lower temperature, this transformation occurs for G of around 
1.8-2.0 e. The simulations which lead to the L-form with the value of 
the torque being close to this limiting value may last for up to 15000r. 
This time becomes significantly shorter for larger values of G. 

The S-form of the DNA was experimentally characterised by over- 
stretching of the dsDNA molecule by about 60%. In models I and II 
the elongation of the system leads to tightening of the bonds between 
p-beads along the chains, which finally leads to significant increase 
in the applied force. Thus the S-form region in the dsDNA phase 
diagram corresponds to structures in which the backbone forms the 
straight line, without imposing the condition of 60% overstretching. 
In T = 0.4 e/A such structures occur while there is applied the force 
of 0.5 e/A. For lower temperature there is needed a bit larger force for 
about 0.05-0.1 e/A. The above values were obtained for applied torque 
G of 1.5 e and 1.7 e respectively for temperatures of 0.4 e and 0.2 e. 

The Pauling form is obtained when both the stretching force and 
the positive twisting torque are large. The smallest value of force 
needed to transform the system into the DNA P-form is 0.25 e/A, 
while the torque applied must be of value around 5.0 e. For larger 
stretching forces, G decreases to 2.1 e at F = 1.5 e/A. In the P-form 
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form, the p-beads come closer together while the remaining beads (b 
and h) become exposed and face out of the helix. 

10 Stretching at constant angular speed 
in scheme D 

In order to study stretching at a constant angular speed, we an- 
chor the bottom beads and attach two frames to the top. Each of these 
frames is as described in the previous section and they coincide ini- 
tially. The beads in one frame are connected to their twins by elastic 
springs. As the outer frame rotates at a constant angular speed, these 
springs get stretched and impose a twist on the inner frame which 
is glued to the DNA. This construct facilitates determination of the 
resistive torque as it is accomplished by monitoring stretching of the 
interframe springs. 

Figure [15] shows the torque of resistance to twisting as a function 
function of the angle of rotation of the outer frame. Two magnitudes of 
the angular speed, u>, are used, 0.00014 1 and 0.00069 ± which differ 
by the factor of 5. We also probe two senses of the twisting: agreeing 
with the helical twist (u > 0) or opposing it (lo < 0). The former 
leads to overtwisting and an indefinite growth in the resistive torque 
due to an increasing infringement of the steric constraints. The latter 
results in unwinding and in a transition from the B-form to the In- 
form. The results clearly depend on the twisting speed. In particular, 
the average torques are 0.319 e and 0.890 e for the smaller and faster 
negative angular speeds respectively. The peaks in the torque result 
from the distortions, but no contacts get ruptured when one infers 
about it from the distance-based criterion. 

11 Conclusions 

The coarse-grained models of the dsDNA discussed here allow for 
studies of features at the level of a single nucleotide. These models 
are found to be fairly equivalent and indicate that sequence-specific 
events can be observed in mechanical manipulations performed at low 
temperatures. However, this capability becomes borderline around the 
room temperature. 

Nevertheless the models proposed here should be useful when 
studying DNA-protein complexes and when assisting nanotechnologi- 
cal DNA assembly processes theoretically. Examples of such processes 
are described in references [401 I4T1 142] . 



15 



Many fruitfull discussions with P. Szymczak are appreciated. This 
work has been supported by the grant N N202 0852 33 from the Min- 
istry of Science and Higher Education in Poland. 



16 



References 

[1] C. Bustamante, Z. Bryant, and S. B. Smith, Nature 421, 423-426 
(2003). 

[2] D. A. Koster, V. Croquette, C. Dekker, S. Shuman, and N. H. 
Dekker, Nature 434, 671-674 (2005). 

[3] J. F. Marko, Proc. Natl. Acad. Sci. (USA) 94, 11770-11772 (1997). 

[4] J. I. Sulkowska and M. Cieplak, J. Phys.: Cond. Mat. 19, 283201 
(2007). 

[5] M. Rief, M. Gautel, F. Oesterhelt, J. M. Fernandez, and H. E. 
Gaub, Science 276 1109-1112 (1997). 

[6] S. B. Fowler, R. B. Best, J. L. Toca Herrera, T. J. Rutherford 
and A. Steward, E. Paci, M. Karplus, and J. Clarke, J. Mol. Biol. 
322, 841-849 (2002). 

[7] G. Yang, C. Cecconi, W. A. Baase, I. R. Vetter, W. A. Breyer, J. 
A. Haack, B. W. Matthews, F. W. Dahlquist, and C. Bustamante 
Proc. Natl. Acad. Sci. (USA) 97, 139-144 (2000). 

[8] M. Carrion- Vazquez, A. F. Oberhauser, T. E. Fisher, P. E. 
Marszalek, H. Li, and J. M. Fernandez, Prog. Biophys. Mol. Biol. 
74, 63-91 (2000). 

[9] H. Lu and K. Schulten, Chemical Physics 247 141-153 (1999). 

[10] E. Paci and M. Karplus, Proc. Natl. Acad. Sci. (USA) 97, 6521- 
6526 (2000). 

[11] G. Pabon and L. M. Amzel, Biophys. J. 91, 467-472 (2006). 

[12] F. Marko and E. D. Siggia, Macromolecules 27, 981 (1994); ibid 
28, 8759 (1995). 

[13] R. M. Jendrejack, J. J. de Pablo, and M. D. Graham, J. Chem. 
Phys. 116 7752-7759 (2002). 

[14] U. Bockelmann, B. Essevaz-Roulet, and F. Heslot, Phys. Rev. 
Lett. 79, 4489-4492 (1997). 

[15] U. Bockelmann, Ph. Thomen, B. Essevaz-Roulet, V. Viasnoff, and 
F. Heslot, Biophys. J. 82, 1537-1553 (2002). 

[16] V. Baldazzi, S. Cocco, E. Marinari, R. Monasson, Phys. Rev. 
Lett. 96, 128102 (2006). 

[17] O. K. Dudko, G. Hummer, and A. Szabo, Proc. Natl. Acad. Sci. 
(USA) 105, 15755-15760 (2008). 

[18] P. Cluzel, A. Lebrun, C. Heller, R. Lavery, J.-L. Viovy, D. Chate- 
nay, and F. Caron, Science 271, 792-794 (1996). 



17 



[19] B. Luan and A. Aksimentiev, Phys. Rev. Lett. 101, 118101 
(2008). 

[20] L. Oroszi, P. Galajda, H. Kirei, S. Bottka, and P. Ormos, Phys. 
Rev. Lett. 97 058301 (2006). 

[21] J. Wereszczynski and I. Andricioaei, Proc. Natl. Acad. Sci. (USA) 
103, 16200-16205 (2006). 

[22] Z. Bryant, M. D. Stone, J. Gore, S. B. Smith, N. R. Cozzarelli, 
and C. Bustamante, Nature 424, 338-341 (2003). 

[23] J. F. Allemand, D. Bensimon, R. Lavery, and V. Croquette, Proc. 
Natl. Acad. Sci. (USA) 95, 14152-14157 (1998). 

[24] B. Albert, D. Bray, J.Lewis, M. Raff, K. Roberts, and J. D. Wat- 
son, Molecular Biology of the Cell, third edition, Garland Pub- 
lishing, New York (1994). 

[25] M. Cieplak and T. X. Hoang, Biophys. J. 84, 475-488 (2003). 

[26] M. Cieplak, T. X. Hoang, and M. O. Robbins, Proteins 49, 114- 
124 (2002). 

[27] M. Cieplak, T. X. Hoang, and M. O. Robbins, Proteins 56, 285- 
297 (2004). 

[28] J. I. Kwiecihska and M. Cieplak, J. Phys.: Cond. Mat. 17, S1565- 
S1580 (2005). 

[29] J. I. Sulkowska and M. Cieplak, Biophys. J. 95, 3174-3191 (2008). 

[30] F. Trovato and V. Tozzini, J. Phys. Chem. B 112, 13197-13200 
(2008). 

[31] T. A. Knotts IV, N. Rathore, D. C. Schwartz, J. J. de Pablo J. 
Chem. Phys. 126, 84901 (2007). 

[32] G. A. Leonard and W. N. Hunter, J. Mol. Biol. 234, 198-208 
(1993). 

[33] C. Clementi, H. Nymeyer, and J. N. Onuchic, J. Mol. Biol. 298, 
937 (2000). 

[34] M. P. Allen and D.J. Tildesley, Computer Simulation of Liquids 
(Oxford University Press, New York, 1987). 

[35] C. Hyeon and D. Thirumalai, Proc. Natl. Acad. Sci. (USA) 102, 
6789-6794 (2005). 

[36] T. E. Ouldridge, I. G. Johnston, A. A. Louis and J. P. K. Doye, 
J. Chem. Phys. (in press) 

[37] N. K. Voulgarakis, A. Redondo, A. R. Bishop, and K. 0. Ras- 
mussen, Phys. Rev. Lett. 96, 248101 (2006). 



18 



[38] M. Cieplak and T. X. Hoang, J. Chem. Phys. 113, 8319-8328 
(2001). 

[39] A. Sarkar, J. F. Leger, D. Chatenay, and J. F. Marko, Phys. Rev. 
E 63, 051903 (2001). 

[40] N. C. Seeman, Nature 421, 427-430 (2003). 

[41] J. C. Mitchell, J. R. Harris, J. Malo, J. Bath, and A. J. Turber- 
field, JACS 126, 16342-16343 (2004). 

[42] R. P. Goodman, I. A. T. Schaap, C. F. Tardin, C. M. Erben, 
R. M. Berry, C. F. Schmidt, and A. J. Turberfield, Science 310, 
1661-1665 (2005). 



19 



Site 


Coordinates\\ 


x y z 


MODEL I 


P 


1 nQO 

— I .Uoy 


— z.zo4 


n a no 
— U.4yz 


h ( A\ 

b {A) 


— o.o41 


— U.blD 


n i in 

U.11U 


U 1 ( A \ 

ni{A ) 


n /i oq 
— U.4yo 


U.ly4 


U.Zob 


hot A \ 

tiZ[A ) 


— u.ozy 


Z.oUo 


1 1 QQ 

l.loo 


b (1 ) 


A A QO 


— U.zyo 


n a sri 
— U.4oy 


hi ( r p\ 


— Z.OZu 


U.014 


n okk 
— U.zoo 


h r >l r r\ 
tiZyl ) 


— l.oyz 




n 71 7 
— U. / 1 1 


6(C) 


-3.737 


-0.997 


-0.186 


61(C) 


-0.793 


2.423 


0.731 


h2(G) 


-0.449 


0.300 


0.131 


63(C) 


-0.005 


-1.851 


-0.500 


6(C) 


-4.603 


-0.802 


-0.163 


hl(C) 


-1.768 


2.928 


-0.103 


h2(C) 


-2.146 


0.665 


0.178 


6.3(C) 


-2.524 


-1.552 


0.428 


MODEL II 


P 


-7.039 


-2.284 


-0.492 


b(A) 


-2.443 


0.261 


0.832 


b(T) 


-3.197 


1.046 


-0.351 


6(C) 


-2.183 


0.230 


0.468 


b(C) 


-3.153 


0.825 


0.139 



Table 1: Cartesian coordinates of the beads in models I and II of the ds- 
DNA. The coordinates depend on the identity of a nucleotide. The follow- 
ing description shows how to generate the DNA double helix for an arbi- 
trary sequence. For the nucleotide, which is placed in the {n + l)-th posi- 
tion in the first strand, coordinates of the related beads can be obtained 
through the following transformation: x(n + 1) = xcos(n 36°) — ysin(n 36°), 
y{n + 1) = xsin{n 36°) + ycos{n 36°), z{n + 1) — z + n3.4, where x, y and 
z denote the starting coordinates as listed in the table. This transformations 
involve rotation around the z axis by 36°, and a shift by 3.4 A along the 
z axis. The second strand in the helix is constructed in a similar way, but 
the initial coordinates have to be transformed from (x,y,z) into (-x,y,-z) for 
every bead belonging to the first nucleotide of the second strand. Then one 
apllies the prior transition to generate the positions of the sites belonging 
to nucleotide paired with {i + l)-th nucleotide in the first strand. A similar 
construction for model III is presented in ref. |31j. 
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Figure 1: Four possibilities of manipulation of the DNA double helix. 
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2 hydrogen bonds 




3 hydrogen bonds 




Figure 2: A schematic representation of model I of the dsDNA. It shows 
formation of 2 and 3 hydrogen bonds. 
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Figure 3: A schematic representation of models II and III. The thin lines 
indicate the way the backbone chains are constructed. 
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DNA atomic structure 



Model I 



Model II 




Figure 4: The atomic representation of the 119D dsDNA structure is shown 
on the left. The remaining panels show the corresponding coarse grained 
representations considered in this paper. 
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Figure 5: The A-type stretching at T = 0. The snapshots at the top show 
subsequent conformations of the pulled dsDNA in model II for u p =0.05 A/r. 
The panels below on the left correspond to v p =0.0b A/r and those on the 
right to f p =0.005 A/r The middle panels show the F — d curves for the three 
models as indicated. At the end of the process, the two chains get fully sep- 
arated and the force drops to 0. The bottom panels show the corresponding 
scenarios of unfolding. 
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Figure 6: Stretching in the A-type mode at T = 0.2e/kB (panels on the left) 
and at T = OAe/ks (the panels on the right). The pulling velocity is 0.05 
A/r. 
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Figure 8: The average force arising during the A-type stretching for three 
different unravelling velocities in model I. The inset shows the log-log plot of 
< F > vs. v n . 
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Figure 9: The B-type stretching for v p =0.05 A/r at two different tempera- 
tures for the three models. The panels on the left represent the results cor- 
responding to T = 0, and the results shown on the right panels correspond 
to T = 0.2e/fce. The snapshots presented at the top show conformations 
during the stretching process of model II at T = 0.0. 
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Figure 10: The C-type stretching for v p =0.05 A/r at two different tempera- 
tures for the three models. The panels on the left represent the results cor- 
responding to T = 0, and the results shown on the right panels correspond 
to T = 0.2e/fce. The snapshots presented at the top show conformations 
during the stretching process of model II at T = 0. 



30 



10 



7.5 



X 
S3 

6 



2.5 


10 



o< 7-5 
1 ' 5 

cd 

G 

fc 2.5 



- 




V = 0.05 A/t " 


— 


\ \\ III 

\ \ \ 


- 




I \ 1 











1 1 


1 1 1 1 


1 1 1 1 1 1 1 1 





0.2 


0.4 0.6 




a\ III 


v = 0.005 A/t ; 




l\ \\ 








"^^s. II 

W A A A ^ " 







J I L 



J I I I l_ 







0.2 0.4 0.6 

T [e/k B ] 



Figure 11: The maximal force for the C-type stretching for v p =0.05 A/t and 
for v p =0.00b A/ t as a function of temperature for the three models. 
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Figure 12: The C-type stretching in model I. The main figure shows F max 
as a function of T for three values of v p . The inset shows the log-log plot of 

Fmax VS . Up . 
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Figure 13: A schematic view of the beads mounted at the DNA end in order 
to apply torque through them. All beads are connected with the springs. 
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Figure 14: The phase diagram representing the final type of dsDNA structure 
obtained after applying the stretching force, F, and the torqe, G, in model 
II. L denotes the L-DNA form, and S signifies stretched chains, in which the 
backbones form straight lines. B denotes the original B-type structure. P 
corresponds to the Pauling form of the DNA, where the backbones get closer 
to each other, and bases stay outside of the helix. The solid lines represent 
results obtained for T = 0.4, while the dashed ones for T = 0.2. 



34 



|w| = 0.00014 1/t 

-cj>0 .,f- f ' 




I I I I I I I I I I L 

250 500 750 



co\ = 0.00069 1/r 




i i I i i I i i I i i 

250 500 750 

angle [degrees] 

Figure 15: Torque of resistance to the twist as a function of the angle of 
rotation for two magnitudes of the angular velocities in T = 0.4 as indicated. 
The system consists of 20 base pairs and the simulations have been performed 
within model II. The dashed lines correspond to a sense of the twist that 
agrees with helical rotation of the dsDNA. The solid lines correspond to the 
opposite sense of the twist. 
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